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Figure  1:  A  HMMWV  traversing  a  deformable  terrain  modeled  using  the  discrete  element  method  (DEM). 
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Abstract 

The  observation  motivating  this  contribution  was  a  perceived  lack  of  expeditious  deformable  terrain  models 
that  can  match  in  mobility  analysis  studies  the  level  of  fidelity  delivered  by  today’s  vehicle  models.  Typi¬ 
cally,  the  deformable  terrain-tire  interaction  has  been  modeled  using  Finite  Element  Method  (FEM),  which 
continues  to  require  prohibitively  long  analysis  times  owing  to  the  complexity  of  soil  behavior.  Recent  at¬ 
tempts  to  model  deformable  terrain  have  resorted  to  the  use  of  the  Discrete  Element  Method  (DEM)  to 
capture  the  soil’s  complex  interaction  with  a  wheeled  vehicle.  We  assess  herein  a  DEM  approach  that  em¬ 
ploys  a  complementarity  condition  to  enforce  non-penetration  between  colliding  rigid  bodies  that  make  up 
the  deformable  terrain.  To  this  end,  we  consider  three  standard  terramechanics  experiments:  direct  shear, 
pressure- sinkage,  and  single-wheel  tests.  We  report  on  the  validation  of  the  complementarity  form  of  con¬ 
tact  dynamics  with  friction,  assess  the  potential  of  the  DEM-based  exploration  of  fundamental  phenomena 
in  terramechanics,  and  identify  numerical  solution  challenges  associated  with  solving  large-scale,  quadratic 
optimization  problems  with  conic  constraints. 
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1.  DEM  for  Terramechanics:  Modeling  and  Numerical  Solution  Strategies  Adopted 

This  contribution  is  motivated  by  an  ongoing  effort  to  identify  predictive  modeling  approaches  that 
can  characterize  dynamics  of  off-road  vehicles.  The  salient  feature  of  off-road  maneuvers  is  the  presence  of 
deformable  terrain.  Owing  to  the  complex  soil  behavior,  deformable  terrain  continues  to  pose  significant 
hurdles  that  limit  the  spectrum  of  scenarios  that  can  be  analyzed  through  computer  modeling.  The  task 
undertaken  is  timely,  given  that  it  is  difficult  and  expensive  to  evaluate  a  vehicle’s  performance  during 
a  majority  of  off-road  maneuvers  using  physical  experiments.  Indeed,  the  range  of  scenarios  that  can  be 
considered  for  physical  testing  is  limited  due  to  time  and  cost  constraints.  It  is  thus  desirable  to  employ 
computer  modeling  in  a  virtual  prototyping  exercise  that,  when  drawing  on  physics-based,  predictive  models, 
can  improve  designs,  compress  the  release  cycle,  and  reduce  costs. 

One  off- road  scenario  of  interest  is  shown  in  Fig.  1.  Therein,  the  vehicle  operates  on  sand,  gravel,  or 
rocky-type  soil,  which  are  modeled  using  DEM.  DEM  represents  soil  as  a  multitude  of  three-dimensional  rigid 
bodies,  called  elements,  where  each  element  is  defined  by  its  size,  shape,  position,  velocity,  and  orientation. 
By  modeling  soil  using  individual  elements,  DEM  allows  for  significant  soil  deformation  and  transport, 
and  the  modification  of  properties  such  as  soil  packing  structure  and  non-homogeneity.  There  are  multiple 
formulations  of  DEM,  classified  based  on  how  the  contact  and  impact  are  handled  when  two  bodies  collide. 
This  paper  assesses  the  predictive  attribute  of  the  complementarity  method,  which  models  contact  as  a 
differential  inclusion.  To  that  end,  it  compares  results  obtained  from  standard  terramechanics  experiments 
to  corresponding  computer  modeling  analyses.  These  experiments  include  direct  shear,  pressure-sinkage, 
and  single  wheel  tests. 

1.1.  Modeling  Frictional  Contact  Via  Differential  Variational  Inequalities 

Consider  a  three  dimensional  (3D)  system  of  rigid  bodies  which  may  interact  through  frictional  contact. 
An  absolute  Cartesian  coordinate  system  will  be  used  to  define  the  equations  of  motion  for  the  time  evolution 
of  such  a  system  [1].  Therefore,  the  generalized  positions  q  =  [rf ,  ef, . . . ,  r^6,  e^J  and  their  time  derivatives 

q  =  •  •  •  ,  r^6,e^JT  are  used  to  describe  the  state  of  the  system.  Here,  Tj  is  the  absolute  position 

of  the  center  of  mass  of  body  j  and  ej  is  the  quaternion  used  to  represent  rotation.  Note  that  the  angular 
velocity  of  body  j  in  local  coordinates,  C #•,  may  be  used  in  place  of  the  time  derivative  of  the  rotation 
quaternion.  Then,  the  vector  of  generalized  velocities  v  =  . . . ,  can  be  related  to  q  via  a 

linear  mapping  given  as  q  —  T  (q)  v  [1]. 

Due  to  the  rigid  body  assumption  and  the  choice  of  centroidal  reference  frames,  the  generalized  mass 
matrix  M  is  constant  and  diagonal.  Further,  let  f  (£,  q,  v)  be  a  set  of  generalized  external  forces  which  act 
on  the  bodies  in  the  system.  Finally,  the  second-order  differential  equations  which  govern  the  time  evolution 
of  the  system  can  be  written  in  the  matrix  form  MV  =  f  (£,  q,  v)  [2]. 

The  rigid  body  assumption  implies  that  elements  that  come  into  contact  should  not  penetrate  each  other. 
Such  a  condition  is  enforced  here  through  unilateral  constraints.  To  enforce  the  non-penetration  constraint, 
a  gap  function,  4>(q,  £),  must  be  defined  for  each  pair  of  near-enough  bodies.  When  two  bodies  are  in 
contact,  or  (q,  t)  —  0,  a  normal  force  acts  on  each  of  the  two  bodies  at  the  contact  point.  When  a  pair  of 
bodies  is  not  in  contact,  or  4>  (q ,£)  >  0,  no  normal  force  exists.  This  captures  a  complementarity  condition, 
where  one  of  two  scenarios  must  hold.  Either  the  gap  is  positive  and  the  normal  force  is  exactly  zero,  or 
vice-versa:  the  gap  is  zero,  and  the  normal  force  is  positive. 

When  a  pair  of  bodies  is  in  contact,  friction  forces  may  be  introduced  into  the  system  through  the 
Coulomb  friction  model  [3].  The  force  associated  with  contact  i  can  then  be  decomposed  into  the  normal 
component,  F ^  —  7z,nn^  and  the  tangential  component,  F ^  +7 where  multipliers  >  0, 

and  7 represent  the  magnitude  of  the  force  in  each  direction,  is  the  normal  direction  at  the  contact 
point,  and  and  w i  are  two  vectors  in  the  contact  plane  such  that  n^,  u^,  and  w i  are  mutually  orthonormal. 

Let  the  contact  points  in  the  local  coordinates  of  each  body  be  expressed  as  and  s respectively.  The 
governing  differential  equations  are  obtained  by  combining  the  rigid  body  dynamics  equations  with  the 
unilateral  constraint  equations.  Then,  the  governing  differential  equations,  which  assume  the  form  of  a 
differential  variational  inequality  (DVI)  problem,  become  [4], 
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q 

=  T  (q)  v 

(1) 

Nc 

M(q)v 

=  f  (*,  q> v)  +  E  hi.n-Dj'n  +  7 i,uDj}U  +  7 i,wDjw) 

A _ 1 

(2) 

0 

2 — 1 

<  (q7)  -L  %n  >  0  %  =  1,2,...,NC 

(3) 

(Ti,it>  li,w) 

=  argmin  (7i,«vTA,u  +  ji,wvTDitW)  . 

(4) 

The  tangent  space  generators  Di  —  [D^n,  Di,u,  Di,w\  €  M6nf>x3  are  defined  as 


Dj  =  [0,  •  -  -  ,  AfpAASi:A,  0,  •  •  •  ,0,  -AjpABsi:Bi  •  *  *  ,  0]  , 


(5) 


and  are  used  to  transform  the  contact  forces  from  local  to  global  frame,  fi{  is  the  coefficient  of  friction 
for  contact  i,  and  Nc  is  the  total  number  of  possible  contacts.  Note  that  for  contact  between  bodies  A 
and  B ,  the  matrix  A^p  —  [n^u^wq]  E  M3x3  is  used  to  represent  the  orientation  of  the  contact  in  global 
coordinates,  and  the  matrices  A  a  and  AB  as  the  rotation  matrices  of  bodies  A  and  £>,  respectively,  where 
the  skew-symmetric  cross-product  matrix  is  defined  as 


s 


-Sz 

0 

$x 


0 


(6) 


1.2.  The  Numerical  Solution  Scheme 

Eq.  1  through  4  are  discretized  to  obtain  an  approximation  of  the  solution  at  discrete  instants  in  time. 
In  the  following,  superscript  (/)  denotes  a  variable  at  time  instant  .  For  example,  q®  and  v®  represent 
the  position  and  velocity  at  time  t^l\  respectively.  Further,  7^  =  h^i  is  the  contact  impulse  for  contact  i. 
Then,  the  discretized  form  of  the  equations  of  motion  is  posed  as  [5]: 


q(m)  = 

M  (v(/+1)  -  v®)  = 


0  < 

hi,u^i,w)  = 


qW  +  hT  (q(0)  v( l+1) 

(7) 

.  iVc 

hi  (*W  qd), 

v(0)  +  ST  (t 'i,nDl„ 
2=1 

i  U  7 i,uD^u  +  r)itWDi  w) 

(8) 

7'(o<0'‘) 

>0,  i  =  1,2,..., iVc 

(9) 

arg  min 

r a  o  . 

(7 iNl+l)T  Di,u 

+  7 Di,W)  > 

(10) 

^2,u+7z%<Mz7z,n 


where  +  h  for  some  integration  time  step  h.  Methods  for  solving  the  above  problem  in  the 

most  general  case  are  lacking  and  various  strategies  are  employed  to  simplify  it.  For  example,  the  friction 
cones  can  be  approximated  by  linear  faceted  pyramids,  leading  to  a  linear  complementarity  problem  (LCP) 
which  can  be  solved  by  pivoting  or  simplex  methods  [6].  However,  these  approaches  in  the  class  of  direct 
methods  can  have  exponential  complexity  in  the  worst-case  [7].  Another  artifact  is  the  lack  of  isotropy, 
given  that  the  friction  cone  is  approximated  by  a  pyramid.  An  alternative  is  to  introduce  a  relaxation  to 
the  complementarity  constraints,  replacing  0  <  i  (q ®,t)  +  _L  7^?n  >  0  with  the  following  [8]: 


°<i^, 


+  DLv 


G+ 1) 


Ti 


(11) 
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As  quantitatively  discussed  in  [9]  and  qualitatively  seen  from  Eq.  11,  the  modification  is  small  when 

[i  or  the  relative  tangential  velocity  at  the  point  of  contact,  \J (.D^v^1))2  +  (E>^v(/+1))2,  are  small. 
Additionally,  it  has  been  shown  in  [5]  that  the  solution  of  the  modified  scheme  will  still  approach  the 
solution  of  the  original  problem  as  the  step-size  tends  to  zero.  An  iterative  method  has  been  developed 
to  solve  this  problem  [10,  11].  It  can  be  shown  that  solving  the  relaxed  discretized  equations  of  motion  is 
equivalent  to  solving  a  Cone  Complementarity  Problem  (CCP)  of  the  form 

Find  for  i  =  1, . . . ,  Nc 

such  that  Ti  3  yf +1^  _L  —  ^ N +  rj  E  Y°  (12) 

where  T *  =  {[x,y,z]T  E  M3 1 \/ y2  +  z2  <  fiix} 
and  Y°  =  {[x,  y ,  z]T  E  M3|x  <  —fiwjy1  +  z2}, 

with  matrix  N  and  vector  r  defined  as 

N  =  DtM1D ,  (13) 

r  =  b  +  DTM-1k,  (14) 

where  bT  =  [bf , . . . ,  b^J  with  hi  =  0,  0]T  E  R3,  and  k  =  Mv ®  +  hi  (t®,  q®,  v®). 

It  can  be  verified  by  considering  the  Karush-Kuhn- Tucker  first-order  necessary  conditions  that  solving 
the  CCP  of  Eq.  12  is  equivalent  to  solving  a  cone-constrained  quadratic  optimization  (CCQO)  problem  [12]. 
For  the  case  with  only  unilateral  constraints  (contacts),  this  optimization  problem  takes  the  form  [12], 

min  q  (7)  =  ^yT TVy  +  rTy  (15) 

subject  to  7^  E  Ti  for  i  =  1,  2, . . . ,  Nc , 

where  7 i  is  the  triplet  of  multipliers  associated  with  contact  i  and  Y i  is  the  friction  cone  of  contact  i.  Note 
also  that  7  =  [7^ ,  7^, . . . ,  7jvc  ] T ■ 

1.3.  Overall  Solution  Methodology  and  Software  Infrastructure 

In  the  adopted  contact  approach,  the  solution  is  produced  at  a  sequence  of  time  steps  t^  <  t ^  <  ...  < 
tW  =  T finai.  Currently,  the  integration  step  size  is  constant;  i.e.,  h  =  ti—ti- 1  is  the  same  for  any  1  <  n  <  N . 
Note  that  a  set  of  initial  conditions  is  provided  at  the  beginning  of  the  analysis;  i.e.,  at  time  t  —  t^°\  The 
solution  is  then  advanced  from  £®  to  as  follows.  The  optimization  problem  posed  in  Eq.  15  is  solved 

to  recover  the  set  of  Lagrange  multipliers  7  that  are  subsequently  used  in  Eq.  8  to  recover  the  new  set  of 
generalized  velocities  Indeed,  note  that  in  the  latter  equation,  provided  7  is  available,  the  velocities 

vd+!)  are  computed  by  multiplying  from  the  left  by  M .  This  matrix  is  constant  and  block  diagonal.  Once 
the  velocity  v^+1)  is  available,  Eq.  7  is  used  to  compute  the  new  position  and  orientation  of  each  part  in 
the  collection  of  bodies  that  make  up  the  system  of  interest.  It  becomes  manifest  that  the  critical  stage  of 
the  solution  methodology  is  solving  the  optimization  of  Eq.  15,  a  process  performed  at  each  time  step  t/.  To 
that  end,  we  draw  on  a  first  order  method  called  the  Accelerated  Projected  Gradient  Descent  (APGD)  [13], 
which  is  briefly  described  in  the  Appendix.  More  than  80%  of  the  entire  simulation  time  is  spent  in  setting 
up  the  optimization  problem  (Eq.  15)  and  solving  it  via  APGD. 

This  numerical  solution  strategy  is  implemented  in  an  open  source  analysis  tool  called  Chrono  [14]. 
Chrono,  which  is  an  multibody  dynamics  engine  released  under  a  BSD3  license,  produces  at  each  time  step 
the  matrix  N  and  vector  r,  an  operation  that  employs  a  collision  detection  task.  For  instance,  in  a  granular 
material  analysis  that  monitors  the  dynamics  of  10000  bodies,  Chrono  performs  at  each  time  step  ti  a 
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collision  detection  that  would  find  approximately  40  000  collision  events.  Given  that  each  contact  event  has 
a  set  of  three  Lagrange  multipliers  associated  with  it,  Chrono  moves  next  to  solve  the  optimization  problem 
of  Eq.  15,  which  calls  for  finding  the  minimum  of  a  cost  function  q(y)  that  depends  on  approximately 
120  000  variables.  Owing  to  the  complementarity-based  formulation  adopted  herein,  the  analysis  proceeds 
at  a  relatively  large  step  size;  i.e. ,  h  «  0.001s.  This  translates  in  performing  1000  collision  detections  and 
solving  1000  optimization  problems  for  each  one  second  of  dynamics  analysis.  To  speed  up  the  solution 
process,  Chrono  resorts  to  parallel  computing  in  each  stage  of  the  analysis  [15]. 

2.  Experimental  Validation 

This  section  reports  on  results  obtained  in  the  process  of  validating  the  solution  methodology  outlined 
in  Section  §1.  To  this  end,  the  Chrono-generated  results  were  compared  against  experimental  data  from 
MIT  [16]  for  three  tests:  direct  shear,  pressure-sinkage,  and  single  wheel.  Quikrete,  a  commercially  available 
concrete  mix,  was  used  in  all  of  the  experiments  for  this  paper.  Quikrete  is  poorly  graded,  having  an  average 
particle  radius  of  approximately  0.4mm  and  bulk  density  pb  =  1.39 gem-3  [17].  Based  on  this  information, 
the  DEM  simulations  used  uniformly-sized  spheres  or  ellipsoids  with  a  major  radius  of  8  mm  (approximately 
20  x  larger  than  the  actual  particle  size)  to  reduce  the  number  of  bodies.  Although  the  exact  material  density 
was  not  measured,  the  Quikrete  particle  density  was  estimated  to  be  pg  =  2.6 gem-3  using  the  following 
equation 


Pb  =  Mg/Vt  =  pgVg/Vt 


based  on  the  volume  of  the  shear  box  enclosure  (V*),  the  total  mass  of  the  granular  material  (Mg),  and  the 
total  volume  of  the  granular  material  ( Vg ). 

2.1.  Direct  Shear  Test 

The  setup  for  the  direct  shear  test  is  shown  in  Fig.  2.  The  test  is  commonly  used  to  measure  the  shear 
strength  properties  of  a  soil,  specifically  the  cohesion,  angle  of  friction,  and  shear  modulus  [18].  A  soil 
sample  is  placed  in  a  shear  box  aligned  under  a  load  cell,  which  applies  a  normal  force  to  the  soil.  The 
top  of  the  shear  box  is  clamped  while  the  lower  half  can  be  moved  in  a  controlled  fashion  by  a  specified 
displacement.  The  horizontal  force  required  to  displace  the  soil  is  measured  to  produce  a  plot  of  the  shear 
stress  as  a  function  of  shear  displacement.  The  shear  box  used  in  the  experiments  had  an  enclosure  that 
contained  granular  material  and  was  approximately  60  x  60  x  60  mm  in  size.  Two  normal  pressures:  16.9 
and  71.4  kPa,  were  tested  for  loosely-packed,  dry  soil.  Four  tests  were  performed  at  each  normal  load  to 


normal  load 


upper  box 


lower  box 
soil  sample 

porous  plate 


Figure  2:  Schematic  of  the  direct  shear  test  (left),  Chrono  simulation  of  the  direct  shear  test  in  the  sheared/final  configuration 
(right). 
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determine  an  experimental  average  and  standard  deviation.  The  3D  DEM  simulations  were  conducted  by 
generating  an  upper  and  lower  shear  box  using  four  box  collision  primitives.  Due  to  the  larger  particle  size 
in  DEM,  a  larger  shear  enclosure  (120  x  120  x  120  mm)  was  used  for  the  simulation.  Approximately  561 
bodies  were  randomly  generated  in  the  interior  of  the  two  shear  boxes  and  allowed  to  settle  under  the  normal 
load  applied  by  an  upper  box  collision  geometry.  After  settling,  the  lower  shear  box  was  translated  at  a 
constant  rate  of  6.6  x  10-4  m  s-1  for  a  distance  of  3  mm  requiring  Tfinai  =4.55  s  for  the  shearing  phase.  Shear 
stress  was  calculated  by  measuring  the  contact  force  on  the  upper  shear  box  and  dividing  by  the  contacted 
shearing  area.  In  general,  the  shear  stress  vs.  displacement  relationship  levels  off  when  the  material  stops 
expanding  or  contracting,  and  when  interparticle  bonds  are  broken.  The  theoretical  state  at  which  the  shear 
stress  remains  constant  while  the  shear  displacement  increases  may  be  called  the  critical  state,  steady  state, 
or  residual  shear  stress.  In  the  case  of  this  paper,  the  residual  shear  stress  is  used  to  study  how  the  behavior 
of  the  direct  shear  simulations  vary  as  a  function  of  the  DEM  parameters,  such  as  the  friction  coefficient  /a. 

2.1.1.  Tuning  of  Solver  Parameters 

Before  calibrating  the  physical  parameters  of  the  model,  a  study  was  performed  to  select  Chrono  solution 
parameters  that  yielded  sufficiently  accurate  results.  Indeed,  if  too  large  of  an  integration  step  size  h  is 
used  or  an  excessively  lax  convergence  stopping  parameter  r  is  adopted,  the  numerical  results  might  not  be 
“converged” ,  thus  compromising  the  predictive  character  of  the  simulation.  To  determine  the  appropriate  h 
and  r,  several  analyses  were  run  using  granular  material  composed  of  uniformly-sized  spheres  with  a  radius 
of  8mm,  a  density  of  2.6 gem-3,  and  a  sliding  friction  coefficient  of  0.5.  The  normal  load  applied  was 
a  =16.9  kPa.  To  determine  the  necessary  step  size  h,  a  tight  solver  tolerance,  r  =  10-2  N,  was  used  while 
running  analyses  at  various  h  values.  The  results  of  this  exercise  are  summarized  in  Fig.  3a.  Conversely, 
given  that  a  step  size  h  —  lx  10-3  s  appeared  sufficiently  small,  a  parametric  sweep  was  carried  out  over  r 
to  gauge  the  sensitivity  of  the  analysis  results  with  respect  to  the  solver  tolerance.  The  outcomes  of  this 
sweep  are  summarized  in  Fig.  3b.  The  results  obtained  elicit  the  following  conclusions:  the  residual  shear 
stress  increases  as  the  solver  tolerance  is  tightened;  the  shape  of  the  numerical  shear  stress  profiles  comes  in 
line  with  expectations;  and  the  results  obtained  are  “converged”  when  h—  lx  10-3  s  and  r  =5xl0-2  N. 


(a)  (b) 

Figure  3:  Residual  shear  stress  for  the  direct  shear  test  with  a  varying  solver  time  step  and  a  fixed  tolerance  of  r  =lxl0-2  N 
(left),  the  residual  shear  stress  for  the  direct  shear  test  with  a  varying  solver  tolerance  and  a  fixed  time  step  of  h  =1  x  10-3  s 
(right). 


2.1.2.  Calibration  of  Model  Parameters 

Since  granular  dynamics  analysis  is  compute  intensive,  one  would  like  to  advance  the  numerical  solution 
at  a  large  step  size  and  take  a  relatively  small  number  of  iterations  to  resolve  each  time  step.  The  outcome 
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of  the  solver  parameter  tuning  stage  described  in  the  previous  section  was  that  the  step  size  and  convergence 
tolerance  required  for  good  results  was  h  —  lx  10-3  s  and  r  =5  x  1CD2  N,  respectively. 

Next,  we  used  these  solver  settings  to  carry  out  a  model  parameter  calibration  aimed  at  determining  the 
friction  coefficient  and  understanding  whether  particle  shape  can  improve  the  quality  of  the  results.  The 
results  for  the  friction  coefficient  study  are  reported  in  Fig.  4a.  The  normal  load  applied  was  a  =16.9  kPa,  the 
granular  material  used  in  the  experiment  had  density  2.6 gem-3  and  the  shape  was  approximately  spherical 
with  radius  8  mm.  With  the  exception  of  fi  —  0.9,  the  numerical  solution  obtained  suggests  that  the  shear 
stress  profile  raises  as  the  sliding  friction  coefficient  increases.  Qualitatively,  the  shape  of  the  numerical 
shear  stress  curves  matches  well  to  experimental  data,  with  a  sliding  friction  coefficient  (i  —  0.5  giving  the 
closest  agreement. 
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Figure  4:  Residual  shear  stress  for  the  direct  shear  test  with  a  varying  friction  coefficient  and  a  constant  shape  ratio  sr  =?  1.0 
(left).  Residual  shear  stress  for  the  direct  shear  test  with  a  varying  shape  ratio  and  a  constant  friction  coefficient  n  —  0.5  (right). 
The  dotted  and  dashed  lines  correspond  to  the  experimental  residual  shear  stress  average  and  standard  deviation,  respectively. 


The  second  model  parameter  tuning  study  was  performed  to  quantify  the  effect  of  particle  shape  on  the 
direct  shear  test  numerical  results.  The  granular  material  was  made  up  of  uniformly-sized  ellipsoids  with  a 
major  radius  of  8  mm  and  a  varying  shape  ratio  from  0.5  -  1.3.  This  ratio  was  defined  as 


rx  =  r,  ry  =  srx  r,  rz  =  r, 
rx  =  r I sr,  ry  =  r,  rz  =  r/sr, 


if  sr  <  1.0 
otherwise 


This  numerical  experiment,  which  was  carried  out  with  /i  —  0.5,  h  =1  x  10-3  s,  r  =5xlO_2N  ,  and 
a  =16.9  kPa,  led  to  the  results  in  Fig.  4b.  Qualitatively,  the  shear  stress  curves  are  not  very  sensitive 
to  the  shape  ratio.  The  shape  of  the  simulated  direct  shear  profiles  matches  quite  well  to  the  experimental 
data. 


2.1.3.  Predictive  Attribute  Assessment 

For  a  normal  load  a  =16.9  kPa,  the  parameter  selection  process  carried  out  in  sections  §2.1.1  and 
§2.1.2  suggests  that  the  numerical  results  come  close  to  experimental  data  when  /x  =  0.5,  h  —  lx  10_3s, 
r  =5x10-2N,  and  sr  =  0.6.  The  predictive  attribute  of  Chrono  is  assessed  next  by  keeping  these  solver 
and  model  parameters  constant  and  modifying  the  experimental  setup;  i.e.,  the  normal  loading  a.  Should 
Chrono  be  predictive,  the  numerical  results  for  a  =71.4 Pa  would  continue  to  be  close  to  corresponding 
experimental  data  measured  for  this  loading  scenario.  Indeed,  the  plot  in  Fig.  5  confirms  that  the  numerical 
results  and  experimental  data  are  similar  for  both  loading  scenarios.  The  direct  shear  simulations  were 
run  on  an  AMD  Opteron  6274  2.2GHz  processor  using  eight  cores.  For  the  normal  load  a  =16.9  kPa,  the 
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(a) 


(b) 


Figure  5:  Shear  stress  vs.  displacement  for  the  direct  shear  test  with  a  =16.9 kPa  (left)  and  a  =71.4 kPa  (right). 


simulation  required  23.2  s  of  compute  time  per  step  (23  200x  slower  than  real  time)  and  used  an  average  of 
4537  APGD  iterations  to  resolve  an  average  of  2112  collisions  per  step.  For  the  normal  load  a  =71.4  kPa, 
the  simulation  required  80.6  s  of  compute  time  per  step  (80  600 x  slower  than  real  time)  and  used  an  average 
of  8925  APGD  iterations  to  resolve  an  average  of  2129  collisions  per  step. 

2.2.  Pressure- Sinkage  Test 

The  pressure-sinkage  test,  shown  in  Fig.  6,  is  used  to  measure  the  load  bearing  properties  of  a  soil. 
To  this  end,  a  plate  is  pushed  with  a  constant  downward  velocity  to  penetrate  a  soil  sample.  Designed  to 
simulate  loading  rates  similar  to  the  ones  exerted  by  a  wheel  during  motion,  the  test  uses  a  load  cell  to 
record  the  force  induced  by  the  plate.  The  soil  bin  used  in  the  experiments  was  150  x  400  x  160  mm  (W 
x  L  x  H)  in  size.  Two  different  plates  were  used  to  penetrate  into  the  loosely-packed,  dry  soil.  Both  of  the 
plates  had  a  length  of  160  mm  and  a  height  of  10  mm.  Each  plate  had  a  different  width;  i.e.,  30  and  50  mm. 
The  plate  penetrated  the  soil  at  10  mm/ s.  Fifteen  experimental  tests  were  performed  at  each  plate  size. 
The  Chrono  tests  were  conducted  by  generating  a  soil  bin  in  which  5377  bodies  were  randomly  generated, 
dropped,  and  allowed  to  settle  under  gravity.  After  settling,  a  plate  was  moved  down  at  a  constant  rate  of 
lOmms-1  for  a  distance  of  30  mm  requiring  Tfinai  =3s  for  the  pressing  phase.  The  pressure  due  to  sinkage 
was  calculated  by  measuring  the  contact  force  on  the  plate  and  dividing  by  the  length  and  width  of  the  plate 
geometry. 

2.2.1.  Tuning  of  Solver  Parameters 

Before  calibrating  the  model  parameters,  a  study  was  performed  to  determine  the  necessary  solver 
step  size  and  termination  criteria  for  the  pressure-sinkage  test.  The  study  was  performed  with  granular 
material  composed  of  uniformly-sized  spheres  with  a  radius  of  8mm,  a  density  of  2.6 gem-3,  and  fi  —  0.5. 
To  determine  the  necessary  time  step,  the  solver  tolerance  was  set  to  r  =5xl0-2  N  and  the  step  size  h 
was  varied  to  determine  its  effect  on  the  numerical  results.  The  outcomes  of  this  study,  carried  out  for 
a  50  mm  plate  width,  are  reported  in  Fig.  7a.  The  numerical  integration  step  size  h  —  lx  10-3  s,  deemed 
appropriate  for  the  shear  test,  yields  good  results  for  the  pressure  sinkage  test  as  well.  Finally,  using  a  step 
size  h—  lx  10-3  s,  a  parametric  sweep  was  carried  out  over  r  to  gauge  the  sensitivity  of  the  analysis  results 
with  respect  to  the  solver  tolerance.  The  outcomes  of  this  sweep  are  summarized  in  Fig.  7b.  Similar  to 
the  analysis  in  §2.1.1,  these  results  suggest  that  a  solver  tolerance  r  =5xl0-2  N  and  a  solver  time  step 
h  —  lx  10-3  s  is  sufficient  to  characterize  the  physics  of  interest  in  the  pressure-sinkage  test. 
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(c) 

Figure  6:  A  photograph  of  the  pressure-sinkage  test  experiment  in  the  initial  configuration  (a),  Chrono  simulation  of  the  pressure- 
sinkage  test  in  the  filled/initial  configuration  (b),  Chrono  simulation  of  the  pressure-sinkage  test  in  the  pressed/final  configuration 
(c).  The  colors  represent  the  relative  magnitude  of  the  linear  velocity  of  the  bodies  (red  =  fast,  blue  =  slow). 


Time  step  [sec] 

(a) 


Figure  7:  Slope  of  the  pressure  vs.  sinkage  curve  for  various  solver  step  sizes  h  (left).  Slope  of  the  pressure  vs.  sinkage  curve 
for  various  solver  tolerances  r  (right). 


2.2.2.  Calibration  of  Model  Parameters 

The  interest  here  is  gauging  the  influence  of  the  friction  coefficient  {i  and  of  the  particle  shape,  controlled 
through  the  coefficient  sr,  on  the  numerical  results  of  the  pressure-sinkage  test.  The  calibration  tests  were 
performed  with  granular  material  composed  of  uniformly-sized  spheres  with  a  radius  of  8  mm,  a  density  of 
2.6 gem-3,  and  a  varying  sliding  friction  coefficient.  In  accordance  with  the  conclusions  reached  in  Section 
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2.2.1,  this  calibration  was  carried  out  using  a  simulation  time  step  h  —  1  x  10  3s  and  a  solver  tolerance 
r  =5x  10-2  N.  The  plate  width  was  50mm. 
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Figure  8:  Slope  of  the  pressure  vs.  sinkage  curve  for  the  pressure-sinkage  test  with  a  varying  friction  coefficient  and  a  constant 
shape  ratio  sr  =  1.0  (left).  Slope  of  the  pressure  vs.  sinkage  curve  for  a  set  of  shape  ratios  sr  and  constant  friction  coefficient  /x  = 
0.5  (right).  The  dotted  and  dashed  lines  correspond  to  the  experimental  residual  shear  stress  average  and  standard  deviation, 
respectively. 


The  results  in  Fig.  8a  suggest  that  the  slope  of  the  pressure-sinkage  curve  obtained  by  Chrono  is  smaller 
than  the  experimental  one  by  roughly  a  factor  of  five.  In  fact,  varying  the  friction  coefficient  has  vary  little 
effect  on  the  slope  of  the  pressure-sinkage  curve.  We  posit  that  this  discrepancy  is  due  to  the  element  dimen¬ 
sion  and/or  the  actual  shape  of  the  quikrete  material  grains.  By  changing  the  shape  of  the  Chrono  elements, 
for  sr  =  0.6  the  experimental  and  simulation  results  match  very  well.  For  fi  —  0.5,  the  sensitivity  of  the 
results  with  respect  to  the  shape  factor  sr  is  summarized  in  Fig.  8b,  for  sr  between  0.5  and  1.3. 

2.2.3.  Predictive  Attribute  Assessment 


Figure  9:  A  comparison  of  the  experimental  and  tuned  pressure-sinkage  profiles  with  a  width  of  50mm(left)  and  a  comparison 
of  the  experimental  and  predicted  pressure-sinkage  profiles  with  a  width  of  30mm(right). 

The  next  set  of  numerical  experiments  used  the  solver  and  model  parameters  identified  in  sections  §2.2.1 
and  §2.2.2  for  a  plate  of  width  50  mm.  Incidentally,  these  parameters  assumed  values  identical  to  the  ones 
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in  the  shear-stress  test.  Specifically,  the  granular  material  was  composed  of  uniformly-sized  spheres  with  a 
radius  of  8  mm,  sr  =0.6,  (i  —  0.5,  r=5xl0-2  N,  and  h  —  lx  10-3  seconds.  Maintaining  these  parameters,  the 
predictive  attribute  of  Chrono  was  assessed  by  comparing  the  numerical  results  to  experimental  data  when 
the  width  of  the  plate  changed  from  30  to  50  mm.  The  results  obtained,  summarized  in  Fig.  9,  confirm  that, 
as  expected,  the  pressure- sinkage  profiles  migrate  to  higher  values  as  the  plate  width  increases.  Moreover,  the 
predicted  values  are  close  to  the  experimental  results  in  Fig.  9b.  The  pressure-sinkage  simulations  were  run 
on  an  AMD  Opteron  6274  2.2GHz  processor  using  eight  cores.  For  a  plate  of  width  50  mm,  the  simulation 
required  27.9  s  of  compute  time  per  simulation  step  (27  900  x  slower  than  real  time)  and  used  an  average  of 
179  APGD  iterations  to  resolve  an  average  of  18  973  collisions  per  step.  For  a  plate  of  width  50  mm,  the 
simulation  required  30.1  s  of  compute  time  per  simulation  step  (30  100  x  slower  than  real  time)  and  used  an 
average  of  199  APGD  iterations  to  resolve  an  average  of  19  058  collisions  per  step. 

2.3.  Single  Wheel  Test  in  Chrono 

A  single  wheel  test  is  used  to  investigate  a  wheel’s  motion  under  controlled  slip  and  normal  loading 
conditions  within  a  confined  soil  bin  of  dimensions  320  x  800  x  150  mm  (W  x  L  x  H).  The  drawbar 
pull,  wheel  torque,  and  sinkage  were  measured  for  a  lug-less  rigid  wheel  for  several  slip  cases  and  loading 
scenarios.  The  wheel  used  in  this  study  had  a  width  of  160  mm  and  a  radius  rw  =130  mm.  To  produce  a 
desired  constant  slip,  the  wheel  was  rolled  on  the  soil  with  a  constant  angular  velocity  of  uo  =0.3 rads-1  and 
a  certain  fixed  translational  velocity  v  based  on  the  slip  defined  as 

v  —  (1.0  —  slip)  uo  rw  . 

The  numerical  tests  in  Chrono  were  conducted  by  simulating  a  soil  bin  in  which  10  790  bodies  were  randomly 
generated  and  allowed  to  settle  under  gravity.  After  settling,  the  wheel  was  rolled  at  the  desired  slip  ratio 
for  T finai  =8  s. 


(a)  (b) 

Figure  10:  A  photograph  of  the  single  wheel  test  experiment  at  MIT’s  Robotic  Mobility  Group  [19]  (left),  Chrono  simulation  of 
the  single  wheel  test  (right). 


Using  the  set  of  solver  and  model  parameters  selected  in  sections  §2.1  and  §2.2,  a  study  was  performed 
to  gauge  whether  Chrono  can  reproduce  the  experimental  results  at  varying  values  of  wheel  slips.  These 
“predictive  attribute”  verification  tests  were  performed  with  granular  material  composed  of  uniformly-sized 
ellipsoids  with  a  major  radius  of  8  mm,  a  density  of  2.6  gem-3,  sr  =0.6,  p  =0.5,  h  =  lx  10-3  seconds, 
r  =5xl0-2  N.  The  quantitative  results  of  this  study  are  summarized  in  Figs.  11  through  13  for  the  drawbar 
pull,  torque,  and  sinkage,  respectively.  It  can  be  seen  that  as  the  slip  of  the  wheel  increases,  the  drawbar 
pull,  torque,  and  sinkage  also  increase  and  the  values  obtained  depend  on  the  normal  loads  applied  to  the 
wheel;  i.e.,  80  N  and  130  N.  Moreover,  the  numerical  values  are  close  to  the  experimental  results.  The  single 
wheel  simulations  were  run  on  an  AMD  Opteron  6274  2.2GHz  processor  using  eight  cores.  For  a  normal 
load  of  80  N,  the  simulation  required  279  s  of  compute  time  per  simulation  step  (279  000 x  slower  than  real 
time)  and  used  an  average  of  473  APGD  iterations  to  resolve  an  average  of  39  090  collisions  per  step.  For 
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a  normal  load  of  130  N,  the  simulation  required  314  s  of  compute  time  per  simulation  step  (314  000 x  slower 
than  real  time)  and  used  an  average  of  513  APGD  iterations  to  resolve  an  average  of  39  379  collisions  per 
step. 


30 


20  |- 
10 
0  - 
-10  - 
-20 
-30 
-40 


-50 


“i - 1 - 1 - 1 - r~ 


-  I 

r  O 


Simulation 
•  Experiment 


-0.8  -0.6  -0.4  -0.2 


0 

Slip 

(a) 


0.2  0.4  0.6  0.8 


60 

i - 1 

1 - 1 

- 1 - 

1 - 1 

1  1 

1 - 

40 

o  - 

20 

I 

2 

I 

Z 

X 

o 

I  ° 

I 

- 

(T3 

-§  -20 

2 

Q 

-40 

I 

I 

o 

O  H 

o 

o 

■ 

-60 

-  i 

o 

o 

" 

-80 

1 - 1 

. 

-0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8 

Slip 

(b) 


Figure  11:  Drawbar  pull  vs.  slip  curves  for  a  normal  load  80  N  (left),  and  130  N  (right). 
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Figure  12:  Torque  vs.  slip  curves  for  a  normal  load  80  N  (left),  and  130  N  (right). 


3.  Conclusions  and  Future  Research 

This  contribution  summarizes  a  validation  study  carried  out  in  order  to  assess  the  predictive  attribute  of 
a  DVI-based  frictional  contact  model  that  is  expected  to  handle  deformable  terrain  problems  in  which  the 
soil  is  represented  using  a  large  number  of  discrete  elements.  The  DVI  model  is  based  on  complementarity 
conditions  for  contact  and  differential  inclusions  for  handling  friction  forces.  The  modeling  methodology  has 
been  implemented  in  an  open  source  dynamics  engine  called  Chrono,  which  is  used  for  soft  soil  ground  vehicle 
mobility  studies.  The  conclusion  of  this  study  is  that  the  predictive  attribute  of  the  modeling  methodology, 
as  exposed  by  its  implementation  in  Chrono,  is  good.  Indeed,  a  unique  set  of  solution  and  model  parameters 
were  used  to  match  experimental  data  in  three  tests:  granular  material  shearing,  pressure  sinkage,  and 
drawbar  pull  at  various  wheel  slip  levels.  The  simulation  times  are  larger  than  hoped  for,  yet  it  is  not  clear 
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Figure  13:  Sinkage  vs.  slip  curves  for  a  normal  load  80  N  (left),  and  130  N  (right). 


whether  a  penalty-based  DEM  approach  would  have  yielded  a  numerical  solution  more  expeditiously.  This 
simulation  performance  aspect,  along  with  a  better  understanding  of  the  sensitivity  of  simulation  results 
with  respect  to  the  shape  and  size  of  the  elements  present  opportunities  for  future  research. 
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Appendices 

A.  Solving  the  Cone-Constrained  Quadratic  Optimization 

In  [20],  Nesterov  developed  a  gradient-descent  method  with  an  improved  convergence  rate  of  O  (l/fc2). 
In  fact,  the  method  in  [20]  was  shown  to  be  an  ‘optimal’  first-order  method  for  smooth  problems  [21]  in 
terms  of  its  performance  among  all  first-order  methods,  up  to  a  constant. 

The  following  set  of  equations  represents  one  iteration  of  the  accelerated  gradient  descent  (AGD)  scheme 
[22].  Note  that  yo  =  xo  G  Mn,  6q  =  1,  q  G  [0, 1]  is  a  tuning  parameter,  and  tk  is  the  step  size  for  the  current 
iteration. 


x/o+i  =yfc-4V/(yfe) 

0k+ 1  solves  9k+1  =  (1  -  9k+ 1)  0l  +  q6k+i 

a.  —  I1  “ 

Pk+1  ~  _i_  n 

Uk  '  Uk+1 

Yk+1  —  xfc+ 1  T  Pk+ 1  (xfc+ 1  —  xfc) 


(16) 

(17) 

(18) 
(19) 


Assume  /  (x)  is  convex  and  Lipschitz  continuous  with  constant  L;  i.e.,  ||V/(x)  —  V/(y)  1 1 2  <  L||x- 
y||2,Vx,y  G  Mn.  Then,  the  method  described  by  Eqs.  16  through  19  converges  for  any  tk  <  1/L.  In  the 
above,  note  that  q  —  1  leads  to  6k  =  1,  fik  —  0,  and  y k  —  x&  for  all  k  >  0,  which  reduces  to  the  gradient 
descent  method.  In  general,  the  parameter  q  can  tune  the  performance  of  the  method  depending  on  the 
specifics  of  the  objective  function  /  (x).  For  example,  if  /  (x)  is  also  strongly  convex,  i.e.,  3/i  >  0  :  /  (x)  > 
/  (x*)  +  (/i/2)  ||x  —  x*| |2,  Vx  G  Mn,  then  the  optimal  value  is  q  —  fi/L ,  which  achieves  a  linear  convergence 
rate.  If  the  objective  function  is  not  strongly  convex,  or  the  strong  convexity  parameter  /i  is  unknown,  then 
it  is  often  assumed  that  q  =  0.  Note  that  the  original  statement  of  the  accelerated  method  in  [20]  had  q  —  0, 
so  the  convergence  rate  of  O  (l //c2)  is  still  valid. 

The  AGD  scheme  can  be  extended  to  constrained  optimization  by  ensuring  that  Eq.  16  does  not  leave 
the  feasible  set.  The  resulting  algorithm,  called  Accelerated  Projected  Gradient  Descent  (APGD)  can  be 
expressed  by  the  following  set  of  computations  to  be  performed  at  each  iteration  k  >  0.  Once  again,  let 
yo  =  xo  G  Mn,  and  9q  —  1. 


Xfc+1  =  n c  (yk  -  tkVf  (yk)) 

(20) 

9k+ 1  solves  6k+1  =  (1  -  9k+ 1)  9k 

(21) 

a  _9k{  1  —  9k ) 

Pk+l  /)2_l  a 

Uk  +  Vk+1 

(22) 

yfc+1  =  Xfc+1  +  /3k- |_1  (x/o+i  -  xfe) 

(23) 

When  /  (x)  is  convex  and  Lipschitz  continuous  with  constant  L,  then  the  method  described  by  Eqs.  20  through 
23  converges  for  any  tk  <  1/L.  An  equivalent  algorithm  was  proved  in  [23]  to  converge  with  the  same 
O  (l/fc2)  rate  as  the  AGD  method.  These  computations  are  performed  until  the  residual,  r,  which  is  defined 
as 


r  =  I  |f  1 12 ,  f  =  — —  (t  —  H/c  (7  —  gd  (AT 7  +  r)))  €  M3Wc,  (24) 

9d 

drops  below  a  specified  tolerance.  This  termination  criteria  is  a  scaled  version  of  the  projected  gradient,  and 
is  designed  to  ensure  that  the  computed  projected  gradient  direction  is  tangent  to  the  constraint  manifold 
at  the  current  iterate.  This  is  desirable  because  it  is  known  that  at  the  optimal  solution,  the  gradient  is 
orthogonal  to  the  constraint  manifold.  Therefore,  it  is  logical  to  use  the  component  of  the  gradient  which 
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is  tangent  to  the  constraint  manifold  as  a  measure  of  error.  To  understand  this  residual,  first  note  that  if 
7  =  7*  is  optimal,  then  II^  (7*  —  Qd  (A/7*  +  r))  =  7*,  so  f  =  0  and  r  —  0  as  expected.  Second,  consider  the 
case  when  7  is  not  optimal.  Then,  it  can  be  verified  that 


n*;  (7  -  gd  (N'y  +  r))  =  7  -  gdf.  (25) 

In  the  preceding,  the  left  hand  side  is  equivalent  to  taking  a  step  of  length  gd  in  the  negative  gradient 
direction  and  projecting  back  to  the  feasible  region.  The  right  hand  side  says  that  the  same  point  can  be 
reached  by  taking  a  step  of  length  gd  in  the  direction  opposite  of  f.  In  the  limit,  as  gd  0,  the  direction  f 
approaches  the  plane  tangent  to  the  constraint  manifold.  Note  that  r  could  be  used  to  measure  convergence 
for  any  value  of  gd,  but  a  small  value  was  used  in  practice  for  the  reasons  just  stated. 


A.l.  APGD  Algorithm 

This  section  states  the  overall  algorithm  obtained  when  applying  the  APGD  method  to  the  problem  in 
Eq.  15.  As  stated,  the  algorithm  includes  an  adaptive  step  size  which  may  both  shrink  and  grow,  an  adaptive 
restart  scheme  based  on  the  gradient,  and  a  fail-back  strategy  to  allow  early  termination  [12]. 

Algorithm  APGD(iV,  r,  r,  Nmax) 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

(10) 
(11) 
(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 
(21) 
(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 


70  =  0nc 
70  =  1  nc 

yo  =  70 

#0  =  1 

r  _  11^(70  -70)112 


Lk 


||70-7o||2 


for  k  !=  0  to  Nmax 
g  Nyk  +  r 
7fc+i  =  n K  (yk  -  tkg ) 

while  7fc+iAr7fc+i+7fc+ir  -  \yINyk+yl*+'gT (7fc+i  -yk)  +  5L*ll7fc+i-y*ll2 
Lk  =  2  Lk 

7  =  7 

7fc+i  =  n ic  (y k  -  tkg) 

endwhile 

ek+1  = 

^k+i  =  6k^tri 

Yk+1  —  lfk+l  +  Pk+l  (7fe+l  —  lk) 
r  =  r  ( 7/c+i) 
if  r  <  cmin 

k‘min  —  k 

7  =  7fc+i 

endif 

if  r  <  r 

break 

endif 

if  gT  (7fc+i  -  7fc)  >  0 
y^+i  7^+1 
0k+l  ~  1 

endif 

Lk  —  0.9L^ 

7  =  7 

endfor 

return  Value  at  time  step  ti+ 1,  7^+1  :=  7  . 
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